Let \(\Omega\) be the unit square \([0,1]^2\). Consider the following problem defined over \(\Omega \times [0,T]\):
\[ \begin{cases} \frac{\partial u}{\partial t}-\Delta u = f &in \ \Omega \times [0,T] \\ \quad u = u_0 &in \ \Omega, \ t=0 \\ \quad u = 0 &on \ \partial \Omega,\ t>0 \end{cases} \]
where \(u_0 = \sin( 2\pi x)\sin(2\pi y)\) is the initial condition, \(f(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y) e^{-t}\) is the forcing term and \(\partial \Omega\) is the boundary of \(\Omega\) where we have prescribed homogeneous Dirichelet boundary condition for all \(t\) grater than \(0\). The exact solution of the previous problem is \(u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y) e^{-t}\). The following window wraps all the steps needed to solve the problem.
## spatio-temporal domain
nodes <- rbind(c(0,0), c(1,0), c(1,1), c(0,1))
edges <- cbind(1:nrow(nodes), c(2:nrow(nodes),1))
domain <- Domain(list(nodes=nodes, edges=edges))
domain <- domain %X% c(0,1)
## build mesh
mesh <- build_mesh(domain, maximum_area = 0.00125, minimum_angle = 20)
# time step
dT = 1e-2
mesh$set_deltaT(dT)
## create Functional Space
Vh <- FunctionSpace(mesh)
## define differential operator in its strong formulation
u <- Function(Vh)
## define differential operator in its strong formulation
Lu <- dt(u) -laplace(u)
## forcing term
f <- function(points,times){
res <- matrix(0, nrow=nrow(points), ncol=length(times))
for( t in 1:length(times)){
res[,t] = (8*pi^2 -1) * sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]) * exp(-times[t])
}
return(res)
}
## initial condition
u0 <- function(points){
return(sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]))
}
## create pde
pde <- Pde(Lu, f)
#> tracemem[0x55bcda60e950 -> 0x55bcd895c298]: Pde Pde eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers withCallingHandlers handle_error process_file <Anonymous> <Anonymous> <Anonymous> <Anonymous> do.call saveRDS withCallingHandlers doTryCatch tryCatchOne tryCatchList doTryCatch tryCatchOne tryCatchList tryCatch
# set initial condition and dirichlet BC
pde$set_initial_condition(u0)
pde$set_boundary_condition(fun= 0.0, type="dirichlet")
## solve problem
pde$solve()
## plot solution
plot(u)